Krylov-implementation of the hybridization expansion impurity solver 
and application to 5-orbital models 
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We present an implementation of the hybridization expansion impurity solver which employs sparse ma- 
trix exact-diagonalization techniques to compute the time evolution of the local Hamiltonian. This method 
avoids computationally expensive matrix-matrix multiplications and becomes advantageous over the conven- 
tional implementation for models with 5 or more orbitals. In particular, this method will allow the systematic 
investigation of 7-orbital systems (lanthanide and actinide compounds) within single-site dynamical mean field 
theory. We illustrate the power and usefulness of our approach with dynamical mean field results for a 5-orbital 
model which captures some aspects of the physics of the iron based superconductors. 

PACS numbers: 02.70.Ss,71.10.Fd,71.30.+h,71.10.Hf 



I. INTRODUCTION 

The development of efficient numerical methods to solve 
quantum impurity problems is an active research area. De- 
mand for powerful and flexible impurity solvers is driven by 
the success of dynamical mean field theory (DMFT), which 
approximates Fermionic lattice problems by self-consistent 
solutions of appropriately defined quantum impurity models. ' 
While impurity models are computationally more tractable 
than lattice models, the desire to include spatial correlations 
via cluster extensions^'-'''* or to treat complicated interaction 
terms in realistic descriptions of multi-orbital systems results 
in considerable computational challenges. 

The multi-site or multi-orbital nature of the most relevant 
impurity models favors Monte Carlo methods. In this area, 
considerable progress has been achieved with the recent de- 
velopment of continuous-time or diagrammatic Monte Carlo 
techniques (CTQMC). The CTQMC algorithms come in two 
flavors. Weak coupling solvers^'*"'^ are based on an expansion 
of the partition function in powers of the interaction terms. 
This is the method of choice for large cluster calculations 
of relatively simple models (such as the one-band Hubbard 
model), because the computational effort scales as the cube 
of the system size. The complementary approach is based 
on an expansion of the partition function in the impurity-bath 
hybridization.^ This so-called hybridization expansion tech- 
nique treats the local interactions exactly and can be applied 
to a wide range of models, including the t-J and Kondo-lattice 
model.^'"' However, since the Hilbert space of the local prob- 
lem grows exponentially with the number of sites or orbitals, 
the computational effort scales exponentially, rather than cu- 
bically with system size. Nevertheless, the flexibility of the 
hybridization-expansion method and the information it can 
provide about the relevant states of the atomic system make 
it a desirable tool in particular for the DMFT study of transi- 
tion metal oxides and actinide compounds. Here, we present 
an implementation of this method which enables the reliable 
simulation of models with up to seven orbitals on present-day 
compute clusters with O(IOO) processors. 

The rest of this paper is organized as follows: Section II 



provides a brief review of the hybridization expansion tech- 
nique in the matrix formulation of Refs. 9,10 and Section III 
discusses the new Krylov-based implementation. We demon- 
strate the accuracy and efficiency of the Krylov approach in 
Section IV, and use it in Section V to compute phase diagrams 
for a "toy model" of the pnictides (a five orbital model with 
almost degenerate bands and relatively large Hund coupling 
term). Section VI is a conclusion and outlook. 



II. HYBRIDIZATION EXPANSION IN THE GENERAL 
MATRIX FORMULATION 

A quantum impurity model describes an atom or molecule 
embedded in some host material with which it can exchange 
electrons. The corresponding Hamiltonian H = H\oc + 
ffmix + i^batii contains three terms; i^ioc = Y.a,i3 '^"'^fPl'^^P + 

J2a J s U°''^-''^'^tpl^tpptp^tps describes the impurity (chemi- 
cal potential, interaction and inter-site/orbital hopping terms), 
Hbath = J2a p ^p^p a^p-a ^ bath of non-intcracting electrons 
whose parameters are fixed by the DMFT self-consistency,' 
and the hybridization term H^^^ = I]a,a',p(^p"'"'^i"p."' + 
h.c.) controls the exchange of electrons between the impurity 
and the bath. Diagrammatic Monte Carlo simulation relies 
on an expansion of the partition function Z = Tr[e^'^^] into 
a series of diagrams and the stochastic sampling of collec- 
tions of these diagrams. For the hybridization expansion,^'^-"' 
we split the Hamiltonian into two parts. Hi = H\oc + i?bath 
and H2 — Hniix, and employ an interaction representation 
in which the time evolution of operators is given by Hi: 
0{t) ~ e^^^Oe^'^^^ . In this interaction representation, the 
partition function can be expressed as a time ordered expo- 
nential, which is then expanded into powers of H2, 
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Equation (1) represents the partition function as a sum 
over all configurations c ~ {ti < ... < r„}, 
n — 0, 1, . . ., Ti £ [0, (3) with weight Wc = 



After the expansion, the time evolution (given by Hi) no 
longer couples the impurity and the bath. It therefore becomes 
possible to integrate out the bath degrees of freedom analyti- 
cally to obtain 



■^bathTrioc 
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The configurations c are now collections of n time argu- 
ments Ti < . . . < r„ corresponding to annihilation oper- 
ators with flavor indices ai, . . . , Q!„ and n time arguments 
t[ < . . . < r/j corresponding to creation operators with fla- 
vor indices a'l, . . . , a^. The element i, j of the matrix AI^^ 
is given by the hybridization function Fa' ,aj (t/ — tj), which 
is defined in terms of the hybridization parameters Vp" " and 
the bath energy levels Given the weights ws, a stochastic 
sampling of all relevant configurations c can be implemented 
using local updates such as the random insertion or removal 
of pairs of creation and annihilation operators. 

For the present purpose, the important thing to note is that 
up to the irrelevant constant Zbath the weights consist of two 
factors: Trioc[. . •] evaluates the imaginary- time evolution of 
the quantum impurity for a given sequence of hybridization 
events, while dct gives the contribution of the bath de- 
grees of freedom which have been integrated out. Using fast 
matrix updates, the determinant ratios for local updates can 
be computed in a time 0{n^). The exponential scaling of the 
algorithm is due to the trace factor. With the exception of 
single-site multi-orbital systems with density-density interac- 
tions (for which the occupation number basis is an eigenbasis 
of Hioc and thus the very efficient segment formulation^ can 
be used), the exponential growth of dim(iJioc) with number 
of sites or orbitals means that the simulation of large systems 
becomes computationally expensive. 

The strategy proposed in Ref. 9 was to evaluate the trace in 
the eigenbasis of the local Hamiltonian. In this basis, the time 
evolution operators e^''^'" become diagonal and can be eval- 
uated easily. On the other hand, the operators ip and which 
are sparse and simple in the occupation number basis, become 
complicated matrices in the eigenbasis of H\oc- To facilitate 
the task of multiplying these operator matrices it is important 
to order the eigenstates according to conserved quantum num- 
bers as explained in Ref. 10. The evaluation of the trace is then 
reduced to block matrix multiplications of the form 



E 



Tr„ 



iO)r 



'(e 



-(T'-r)H|„. 



(e— ) 



(3) 

where O is either a creation or annihilation operator, m de- 
notes the index of the matrix block, and the sum runs over 
those sectors which are compatible with the operator se- 
quence. With this technique, 3-orbital models or 4-site clus- 
ters can be simulated efficiently. ' 1.12,13,14,15 However, since the 



matrix blocks are dense and the largest blocks grow exponen- 
tially with system size, the simulation of 5-orbital models be- 
comes already quite expensive and the simulation of 7-orbital 
models with 5, 6 or 7 electrons is only doable if the size of the 
blocks is severely truncated. 

In fact, one should distinguish two types of truncations: 

(i) the truncation of the outer trace (X]^o„t, ,„) to those quan- 
tum number sectors or states which give the dominant 
contribution, 

(ii) the reduction of the size of the operator blocks (0),„' , m" 
via elimination of high-energy states. 

The truncation of type (i) is harmless at low enough temper- 
ature, because it restricts the possible states at only a single 
point on the imaginary-time interval. On the other hand, trun- 
cations of the type (ii), if not done properly, can lead to sys- 
tematic errors, whose effect will be hard to estimate in large 
systems, because the truncations are necessarily severe. 



III. KRYLOV-SPACE METHOD 

As an alternative strategy to evaluate the trace factor in 
Eq. (2) we propose to 

1. adopt the occupation number basis, in which the ip- 
operator matrices can easily be applied to any given 
state, and in which the sparse nature of H\oc can be ex- 
ploited during the imaginary time evolutions by relying 
on efficient Krylov-space methods, 

2. to approximate the outer trace by a sum over the lowest 
energy states (i.e. truncation type (i) introduced above). 

This implementation involves only matrix-vector multiplica- 
tions of the type ip^'^^v) and H\oc\v), with sparse operators 
and i^ioc, and is thus doable in principle even for systems 
for which the multiplication of dense matrix blocks becomes 
prohibitively expensive, or for which the matrix blocks will 
not even fit into the memory anymore. Furthermore, no ap- 
proximation of type (ii) is required, so that all excited states 
remain accessible at intermediate t in the trace. The sparse 
nature of the hybridization operators is evident given the fact 
that they consist of creation and annihilation operators in the 
occupation number basis, ifioc is sparse because the number 
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of interaction terms is proportional to a small integer power 
of the number of orbitals, while the dimension of the matrix 
grows exponentially with the number of orbitals. 

Our implementation is based on very efficient sparse 
matrix algorithms for the evaluation of matrix expo- 
nentials applied to a vector, i.e. exp(— Ti7ioc)|w).'^''^''^ 
These algorithms construct the Krylov space ICp{\v)) ~ 
span{|u), i?ioc|w), i^ioJw), . • • , H^^^lv)} and then approxi- 
mate the full matrix exponential by the matrix exponential 
of the Hamiltonian projected onto the Krylov space ICp{\v)). 
In Ref. 17 it has been shown rigorously that these Krylov 
space algorithms converge rapidly as a function of p, typi- 
cally reaching convergence for very small iteration numbers 
p ^ -/Vdiin, although the number of iterations depends on the 
time interval r. 

Let us describe the algorithm for the trace evaluation in 
some more detail. First, during the initialization part of the 
simulation, the following steps are required: 

1 . Obtain the low energy spectrum and eigenfunctions of 
H\oc using (Band-)Lanczos or Davidson techniques, or 
alternatively diagonalize iJioc completely using full di- 
agonalization techniques. The Band-Lanczos or David- 
son techniques are needed to resolve the exact degen- 
eracies of the eigenfunctions. 

2. Decide which eigenstates of the spectrum are to be kept 
in the outer trace. It is important not to destroy the mul- 
tiplet structure of iJioc when truncating the trace. The 
truncation criteria employed in our implementation are 
discussed in more detail in Sec. IV. 

Then, in the actual evaluation of a trace, we proceed as fol- 
lows: 

3. Propagate a retained state in the trace up to the first hy- 
bridization event (forward and backward in time). Since 
the initial state is an eigenstate of i^ioc, this state is sim- 
ply multiplied by an exponential factor for the first in- 
terval. 

4. Apply the hybridization operator on the propagated 
state. 

5. Propagate the current state up to the next hybridization 
event using the Krylov-space approach to the matrix ex- 
ponential described above. The state to be propagated 
is generically not an eigenstate of i^ioc anymore, so the 
Krylov space must be constructed up to a certain dimen- 
sion. The Krylov space size should not be kept fixed, 
but should be determined for each imaginary time inter- 
val according to a convergence criterion. In the applica- 
tions reported in the present paper the average Krylov 
space dimension is « 2. 

6. Go back to step 4 if more hybridization operators are 
present. 

7. Add the contribution of the propagated state to the trace. 

8. Go back to step 3 until all retained states have been con- 
sidered in the trace. 



In the truncated trace approach it is important to measure 
the various local observables at r = /3/2 in order to be least 
affected by the truncation of the trace at r = (and equiva- 
lently at t = /?). 

We conclude this section by illustrating the main advan- 
tage of the Krylov space method through a simple time com- 
plexity analysis of the algorithm. Say we want to determine 
the trace of a given sequence of the hybridization operators i/; 
and According to the truncation (i) introduced above we 
perform the trace over TVtr < A^dim states, where N^im is the 
typical size of the impurity Hilbert space, which grows ex- 
ponentially with the number of sites or orbitals contained in 
the "impurity". Since there are A^hyb hybridization events, the 
complexity of the application of the hybridization operators 
is 0(A^hyb X A'dim X A'tr). The imaginary time evolution on 
the other hand is nontrivial on A'intervai = Ahyb — 1 intervals. 
Based on Ref. 17, we assume a typical number of iterations 
A'iters <C Adiin is needed to reach convergence for the imagi- 
nary time evolution of a single state \ v) over an interval length 
T. It follows that the complexity of the imaginary time evolu- 
tion part is 0{Nty x A'imeivai x iViter X N^im) and the overall 
time complexity amounts to 

0{Niin, X Nn- X [Nhyb + A^interval X A^i.e,-]). 

In the worst case where we retain all states in the trace Na- = 
Adim the complexity scales as N^^^, but in the best case Atr = 
0(1) the time complexity is linear in Adim- 

In comparison the matrix formulation has a less favorable 
scaling with Adim- In the case where we keep all states in the 
trace the time complexity is O(A'intervaiA'dim) because of the 
expensive dense matrix-matrix multiplications, whereas it is 
0(A'intervai X N^i^ X A't, ) for the truncated trace version. 

While it is therefore obvious that in theory the Krylov space 
approach is the method of choice due to its superior A'dim scal- 
ing, in practice the precise numbers of A^tr, A'iter, and Adim will 
determine which one of the two formulations performs better 
for a given problem with tractable Hilbert space size. In the 
following section we address the performance and scaling of 
the two algorithmic formulations. 

IV. EFFECT OF TRUNCATION AND EFFICIENCY 

A. Accuracy of the Krylov approach 

Before trying to determine the system size for which the 
Krylov implementation outperforms the matrix method, we 
demonstrate the accuracy of the new approach. We consider 
multi-orbital models with a local Hamiltonian of the form 

a, a a 

- X! '^(^li^lt'^^-LV'a,! + V'l,|V'l,xV'a,TV'aa + h.C.) 

(4) 
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FIG. 1 ; (color online) Comparison between the Green functions of a 
3-orbital model computed with the matrix method (open symbols, no 
truncation of the trace) and the Krylov method (full symbols, trun- 
cation of the trace to the lowest energy states) for different temper- 
atures (top panel, f/ = 6) and interaction strengths (bottom panel, 
/? — 50). The results become indistinguishable at temperatures 
which are < 1% of the bandwidth. 
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FIG. 2: (color online) Efficiency (number of updates per second) of 
the different implementations as a function of system size. The mod- 
els are n-orbital impurity models with rotationally invariant Hund 
coupling at half-filling (jj, = /ihaif), U — 6, J/U = 1 /6, (3 — 50. For 
the matrix method we show results without truncation (diamonds) 
and with truncation of the trace to the quantum number sector con- 
taining the ground state (triangles). In the Krylov calculation, the 
trace is truncated to the lowest energy states. 



case, since the ground state of the half-filled model carries 
spin 3/2), which means that we use here the 0(1) approxima- 
tion for iVti. While deviations between the approximate and 
exact result are apparent at (3 = 3.125, they become smaller 
as temperature is lowered and for /3 > 50 can be considered 
negligible. The bottom panel shows that for f3 = 50, essen- 
tially perfect agreement between the two methods is found for 
all relevant interaction strengths. 



and rotationally invariant interactions ([/' = U — 2 J, so the 
inter-orbital interactions are U — 2 J for opposite spin and U — 
3 J for same spin). Figure 1 compares the Green functions for 
a 3-orbital model with Hund coupling parameter J = U/6. 
The hybridization function is that of a non-interacting model 
with semi-circular density of states of bandwidth 4 eV, and the 
chemical potential has been chosen such that the system is at 
half-filling (fi = |J7 — 5 J). The crystal field splittings Aq are 
zero. We give the parameters U, J, fi and A in units of eV. 

For three orbitals, both methods yield accurate results in a 
few CPU hours. The top panel of Fig. (1) shows the mea- 
sured Green functions for J7 = 6, J = C//6 = 1 and differ- 
ent values of inverse temperature /3, while the bottom panel 
compares the results for /3 = 50 and different values of the in- 
teraction strength. The open symbols (red online) were com- 
puted with the matrix method without any truncations. These 
are thus exact results (Monte Carlo errors are much smaller 
than the symbol size) which may be used to test the accuracy 
of the Krylov approach. The Krylov results are plotted by 
full symbols (blue online). These data were computed with 
only the lowest energy states in the outer trace (four in this 



These results can readily be understood from the scaling 
of the perturbation order with U and /? in the hybridization 
expansion method.^ The average perturbation order grows 
roughly linearly with the length /3 of the imaginary time inter- 
val and decreases as interaction strength is increased. Since 
we restrict the system to the ground state at one point of the 
imaginary time interval (t = 0), a larger number of hybridiza- 
tion events facilitates the relaxation into the true equilibrium 
distribution (measurements are performed at r = /3/2). More 
importantly, the lower the temperature, the larger the overlap 
of this probability distribution with the ground state, i.e. the 
probability of the system being in the ground state at any given 
time becomes large. Thus, forcing the system into the ground 
state at r = to compute the trace more efficiently has no se- 
vere effects at low enough temperature. Our data suggest that 
the truncation of to the ground state vectors is legitimate 
for temperatures which are < 1% of the bandwidth (4 eV) and 
we will use this 0(1) truncation in all subsequent Krylov cal- 
culations. We also note that the truncation of the trace does 
not seem to induce a sign problem for the multi-orbital prob- 
lems studied in the present work. 
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FIG. 3: (color online) Efficiency as a function of chemical potential 
for the 5 orbital model (U = 6, J/U = 1/6, /3 = 50). The chem- 
ical potentials have been chosen such that the ground state of H\oc 
has 5, 6, 7, 8, 9 and 10 electrons, respectively. The corresponding 
degeneracy A'tr is plotted next to the Krylov data. 
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FIG. 4: (color online) Weight of the different (n| , rij ) quantum num- 
ber sectors for n — ^haif = 4.98 (M, — 6), fj, — /ihaif = 5 (A^t,- = 31) 
and fi - ^haif = 5.02 (iV„ = 25). 



B. Efficiency 

To compare the efficiency of the two implementations we 
plot in Fig. 2 the number of local updates per second for multi- 
orbital systems with n = 2, 3, . . . orbitals. A local update is 
either an insertion or a removal of a pair of ip, operators 
and involves the calculation of Trioc- In our n-orbital mod- 
els [Eq. (4)] each orbital interacts with every other through 
density-density, spin exchange and pair hopping terms. The 
intra-orbital repulsion is U, the Hund coupling parameter 
J = U/6, and the crystal field splittings are zero. We chose 
L/ = 6, /3 = 50 in all the calculations, and the hybridization 
function of the non-interacting model with semi-circular den- 
sity of states of bandwidth 4. The half-filling condition for 
these multi-orbital systems is /ihaif = — — (n — 1)| J. 
The blue lines with diamonds show the results for the Matrix 
code without any truncation. For n > 3, the evaluation of 
the trace becomes the bottleneck of the simulation and we ob- 
serve an exponential decrease in the number of updates per 
second. The red lines with triangles show the result for the 
matrix code in which the trace is restricted to the sector m 
containing the ground state (but without any truncation in the 
size of the blocks). The rather modest effect of the truncation 
is due to the fact that at half-filling the largest blocks cannot 
be discarded. 

The black lines with circles show the number of updates 
obtained with the Krylov-method. The curve still drops ex- 
ponentially with increasing n, but the slope is smaller than in 
the matrix case, as expected from the scaling argument in the 
previous section. While the number of updates in the matrix 
implementation drops by about 4 orders of magnitude as n is 
increased from 3 to 5, it drops only 2 orders of magnitude in 
the Krylov implementation. The more favorable scaling in the 
Krylov-case allows us to measure also n ~ 6 and n = 7 and as 
seen in Fig. 2, the slope remains essentially unchanged. The 
time per update increases by about two orders of magnitude 



from 71 = 5 to n = 7. Given that the Krylov code allows the 
simulation of 5-orbital models (transition metal compounds) 
on a small number of processors, we therefore expect that 
this method will enable the controlled and accurate simula- 
tion of 7-orbital models (lanthanide and actinide compounds) 
on larger clusters with a few hundred processors. 



C. Effect of the ground state degeneracy 

The outer trace must at least contain all the ground state 
eigenvectors, and the ground state degeneracy of iJioc de- 
pends on the model parameters. We therefore compare in 
Fig. 3 the efficiency of the Matrix and Krylov implemen- 
tations for the five orbital model at different values of the 
chemical potential (chosen such that the ground state lies in 
the ritot = 5, 6, . . . , 10 electron sector). The corresponding 
ground state degeneracies are 6, 25, 40, 30, 10 and 1. The 
increase in N^j from 6 to 25 leads to a slight decrease in the 
efficiency of the Krylov method if jj. is increased from /ihaif to 
A'haif + 6. For even larger /i, the efficiency increases, because 
the relevant quantum sectors become smaller. This also ex- 
plains why the advantage of the Krylov implementation over 
the Matrix implementation decreases as one moves away from 
half-filling. 

If the Krylov trace is restricted to ground state vectors, level 
crossings in H\oc will typically lead to sudden changes in the 
number and types of states considered in Trioc, even in situa- 
tions where the physical state of the system is not expected to 
change dramatically. However, as shown in Fig. 4, this does 
not lead to inconsistencies if the temperature is sufficiently 
low. The figure plots the probability distribution of the dif- 
ferent quantum number sectors, ordered as (7i| = 0, tt.^ = 
0), (n^ = l,ni = 0), . . . , (n^ = 4, ?i| = 5), (t^^ = 5, tt.^ = 
5) from left to right, for the 5 orbital model with [/ = 6, 
J = U/6 and ^ = /^haif + 4.98 (red circles, = 6), 

— Mhaif + 5 (blue diamonds, A^tr = 31) and fi = /Lthaif + 5.02 
(black triangles, N^- — 25) at jS — 50. All three trace calcula- 
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tions yield consistent distributions, and thus the same physical 
state. The (small) inaccuracies near level crossings could be 
further reduced by retaining all the states in a certain energy 
window above the ground state. 



V. APPLICATION 

In this section we illustrate the usefulness and efficiency of 
the Krylov method with DMFT results for 5-orbital models 
with semi-circular density of states of bandwidth 4. We will 
consider the situation in which all bands are degenerate ( Aa = 
0, a = 1, . . . , 5) and a "2+3" eg-t2g crystal field splitting of 
magnitude 0.5, in which the doublet is shifted down (Ai = 
A2 = 0.5, A3 = A4 = A5 = 0). All the calculations are for 
/? = 50 and require less than 50 CPU hours per iteration. 



A. Orbitally degenerate case 

Figure 5 shows the paramagnetic phase diagram in the 
space of chemical potential (relative to ^haif = 4.5f7 — 10 J) 
and interaction strength. The top panel is the result for 
J/U =1/4 and the bottom panel for J/U = 1/6. We first 
discuss the orbitally symmetric case which corresponds to the 
blue lines with stars. The figure shows the Mott insulating 
lobes with n = 5 and 6 electrons. Additional lobes with 
n = 7, ... ,9 and a band insulating solution with n = 10 
also exist, but are not shown (computations near half-filling 
are the most challenging ones, because they involve quantum 
number sectors with high dimension). The Hund coupling 
J is seen to stabilize the half-filled 71 = 5 Mott lobe while 
pushing the critical interaction for the n = 6 Mott lobe to 
larger values. The different widths of these lobes and their J- 
dependence can be understood from the /i-dependence of the 
eigenstates of Hio^ as explained in the context of a 3 -orbital 
calculation in Ref. 12. A comparison of the 5-orbital result 
for J/U ~ 1/6 to the lower panel of Fig. 2 in Ref. 12 and the 
2-orbital calculations (Fig. 2) in Ref. 19 furthermore shows 
the evolution of Uc with increasing number of orbitals; in 
the 2-orbital model f/haiffiUed ^ 3 7^ ^^le 3-orbital model 

jjhalf filled ^ 3 ^jjj [/halffiUed+l ^ jjj (j^g S.Qrbital 

case we find [/half fiUed ' ^ 2 and [/half fiiied+i ^ j^ius, in 
the presence of a Hund coupling J = U/6, the critical inter- 
action strength for the Mott insulating phase decreases with 
increasing number of orbitals, in contrast to the situation for 
J ^ 0.20.21 

The very large interaction strengths required to study Mott 
physics away from half-filling are not a problem. An attrac- 
tive feature of the hybridization expansion method is the fact 
that the relevant perturbation orders decrease with increasing 
interaction strength.^ At [/ = 16 and /3 = 50, the n = 6 Mott 
insulating solution for J/U ~ 1/6 has average perturbation 
order « 1.8 per orbital and spin, i.e. 18 in total. 
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FIG. 5: (color online) Phasediagram of the 5-orbital model in the 
space of chemical potential and interaction strength for J/U — 1/4 
(top panel) and J/U = 1/6 (bottom panel). The blue lines with 
stars show the n — 5 and n — 6 Mott lobes for the model without 
crystal field splitting. The red line with crosses shows the effect of a 
2+3 crystal field splitting of magnitude A = 0.5 (2 orbitals shifted 
down) on the Mott lobes. Both Mott lobes are now contained in an 
orbital selective Mott phase (boundary marked with black circles) in 
which the three degenerate bands are insulating and half-filled, while 
the two degenerate bands are metallic. Error bars are on the order of 
the symbol size. 



B. Crystal field splitting and orbital selective Mott transition 

We now consider the effect of shifting two orbitals down by 
Ai = A2 = 0.5. The resulting phase diagram is plotted with 
crosses and red lines. The n — 5 and 6 Mott lobes are little 
affected by the crystal field spUtting. The value of Uc is almost 
unchanged for n = 5, and decreases by about 2 for n = 6. The 
width of the n = 6 lobe is increased by about A and shrinks 
by a similar amount for n = 5. Compared to the 2+1 splitting 
considered in Ref. 12 the stability of the n — 6 lobe is not 
dramatically enhanced, because the insulating phase does not 
consist of half-filled and band-insulating solutions: the two 
degenerate bands (a = 1, 2) accommodate 3 electrons. Larger 
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FIG. 6: (color online) Expectation value of as a function of J/U 
for the five orbital model with n = 6, U — 2 and (3 — 50. Circles 
show the result for degenerate orbitals, diamonds for a crystal field 
splitting A = 0.5 which shifts two orbitals down. 



effects on the stability of the n ~ 6 lobe are expected for 4+1 
or 3+1+1 splittings. 

A qualitative difference to the orbitally symmetric case is 
that the n = 5 and 6 lobes are now embedded in an orbital 
selective Mott phase characterized by insulating, half-filled 
Green functions in the three degenerate bands (a = 3, 4, 5) 
and metallic Green functions in the two lower lying bands 
(a ~ 1,2). For n ~ 6 electrons, the transition from the 
metallic into the orbital selective Mott phase takes place near 
U 12 (J/U = 1/4) and U « 4.2 (J/U = 1/6), respec- 
tively. The metal-insulator transition in the three degenerate 
bands thus occurs at an interaction strength which is substan- 
tially smaller than the Uc required to induce the Mott transi- 
tion in the model without crystal field splitting. This finding 
is consistent with the enhanced stability of the half-filled Mott 
lobe in the 3-orbital model. 



C. Total spin 

The Krylov implementation retains the attractive features 
of the hybridization expansion, such as the ability to measure 
the relevance of different states in the Hilbert space of H\oc 
(Fig. 4). As a practical application we plot in Fig. 6 the ex- 
pectation value of the total spin squared, 5^, as a function of 
Hund coupling J/U in the weakly correlated metallic phase 
(U ~ 2, (3 ^ 50) with 6 electrons. Results for degenerate 
bands and for a 2+3 crystal field splitting A = 0.5 are shown. 
The atomic ground states for A = (0.5) correspond to 
52 = 0, 2, 6 (0, 2) for J/U = 0, 5^ = 6 (2) for J/U = 0.05 
and = 6 (6) for J/U > 0.1. The atomic picture is how- 
ever not a good reference in the parameter regime considered 
here. We conclude from Fig. 6 that in the moderately cor- 
related metallic phase the lower spin states have appreciable 
weight, and the effect of the crystal field splitting is small. No 



dramatic increase in S*^ is observed as J/U is increased from 
to 0.05, and the crystal field splitting of A = 0.5 leads to no 
significant reduction in S"^, even at J = 0. 



D. Implications for pnictides 

The toy model considered here captures some aspects of the 
iron-based high temperature superconductors. A minimal de- 
scription of these materials seems to require all five d-bands,^^ 
and the bandwidth of 4 eV adopted here is consistent with the 
band structure obtained from density functional theory. Crys- 
tal field splittings appear to be small (A w 0.2-0.5), although 
(due to the tetrahedral coordination) not of the simple 2+3 
type considered in Fig. Sp'^^ No consensus has yet emerged 
about the interaction parameters U and J and the role of cor- 
relations. Some authors argue that U should be quite large 
and the material close to a Mott transition^^ or to an orbitally 
selective Mott state. ^"^'^^ Other theoretical studies, however, 
adopt small interaction parameters, U ^ 2, and J = 0.2- 

0.6.22-26 

Our phase diagrams for J/U ~ 1/6 and J/U = 1/4 
in Fig. 5 can provide some information about the role of 
correlations and crystal field splittings, and the relevance of 
Mott physics. In particular, we see that the solution with 
71 = 6 and U = 2 is far away from the Mott lobe in the 
orbitally symmetric case, especially if the Hund coupling is 
large. On the other hand, our results for the 2+3 splitting 
indicate that even relatively small crystal field splittings can 
have a substantial effect on the phase diagram and lead to the 
opening of a gap in some bands at U much below the Uc for 
the fully gapped phase. Correlations in pnictides may thus be 
relevant in the sense that the materials are not too far from 
an orbital selective Mott state. But for such a scenario, the 
Hund coupling parameter J would have to be rather small: 
for J/U = 1/4:, U = 2 is a factor of 6 below the critical value 
for the orbital selective transition and for J/U = 1/6 it is still 
a factor of 2 below. A moderately correlated metallic state 
seems more consistent with the phase diagram of our simple 
toy model. Figure 6, on the other hand, indicates that the 
total spin in the realistic parameter regime is not particularly 
sensitive to the values of J and A. At [/ = 2, the spin 2 states 
start to dominate only for J > 0.6. However, since details 
of the band structure appear to affect the properties of iron 
based superconductors in profound ways,26 more realistic 
calculations within the LDA+DMFT framework would be 
required to settle these issues. 



VI. CONCLUSIONS 

We have presented an implementation of the hybridization 
expansion impurity solver which makes use of sparse-matrix 
exact diagonalization techniques to evaluate the weight of 
Monte Carlo configurations. This method, while still scal- 
ing exponentially with system size, enables a much more ef- 
ficient simulation of large multi-orbital problems than the es- 
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tablished matrix formulation. In the new approach, the trace is 
restricted to a small number of states (typically the lowest en- 
ergy eigenstates of the Hamiltonian) but no other truncations 
or approximations are necessary during the time evolution of 
these states. We have demonstrated that the restriction of the 
outer trace leads to negligible systematic errors at low tem- 
perature. Therefore, the Krylov method provides a controlled 
and efficient implementation of the hybridization expansion 
approach which enables the DMFT study of transition metal 
and actinide compounds with realistic interactions. 
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